%clc;
%% 生成激励信号
    % ft = 1.0;
    Fs = 48000;  % 采样率：1 GHz
    freqs  = [18000 20000 21000 22000 23000 24000];  % 混合频率
    amps   = [1.00 1.00 1.00 1.00 1.00 1.00];
    pn     = 32768  ;
    % ==== 时间轴 ====
    t = 0 : 1/Fs : ((pn-1)/Fs);
    % ==== 合成信号 ====
    sin_wave = zeros(1,pn);
    for i = 1:length(freqs)
        sin_wave = sin_wave + amps(i) * sin(2 * pi * freqs(i) * t );
    end
    M=length(sin_wave);
    fft_plot_tittle=sprintf("原始数据");
    plotCustomSpectrum(sin_wave(pn/4+1:pn),pn*3/4,Fs,fft_plot_tittle); %% 对原始数据进行FFT之前进行了截取,以避免频谱泄露
%% OSR 设置
    osr = 128;
    fft_osr = 128;

%% Kaiser窗FIR滤波器设计
    fir_osr =  1 ;
    fc = 20000;             % 截止频率20kHz
    tband = 80000;          % 过渡带宽50kHz
    beta = 14.26526;        % Kaiser窗形状参数
    N=200;                   % 滤波器阶数
    % N = round(10*Fs*osr/(2*tband));



%% 一段 fir 2倍过采样插值/滤波
    % 设计一段FIR滤波器参数
    fir_h   = fir1(N, fc/(Fs*fir_osr/2), 'low', kaiser(N+1, beta));

    % % 分析频率响应
    figure();
    freqz(fir_h, 1, 2048, Fs*fir_osr);
    title(['Kaiser窗低通FIR: ', num2str(fc/1000), 'kHz截止, N=', num2str(N)]);

    fft_osr = fir_osr;

    %插值
    fir_inp_data = zeros(1, length(sin_wave)*fir_osr);  % 插值后总长度
    for i=0 : (length(sin_wave)-1)
        fir_inp_data(1+i*fir_osr) = sin_wave(i+1);
    end
    % fft_plot_tittle=sprintf("一段fir2倍插值滤波输入数据");
    % plotCustomSpectrum(fir_inp_data(pn*fir_osr/4+1 : pn*fir_osr),length(fir_inp_data)*3/4,Fs*fir_osr,fft_plot_tittle);%% 对数据进行FFT之前进行了截取,以避免频谱泄露

    % FIR滤波
    fir_1th_y=filter(fir_h,1,fir_inp_data);
    fir_1th_y=fir_1th_y*(fir_osr);
    fft_plot_tittle=sprintf("一段fir2倍插值滤波输出数据,OSR=%d\n",fft_osr);
    plotCustomSpectrum(fir_1th_y(pn*fft_osr/4+1 : pn*fft_osr),length(fir_1th_y)*3/4,Fs*fft_osr,fft_plot_tittle);%% 对数据进行FFT之前进行了截取,以避免频谱泄露

%% 线性插值滤波器设计
    line_input = fir_1th_y;
    line_osr      = osr / (fir_osr); % 过采样率
    fft_osr = fir_osr * line_osr;

    line_y = zeros(1,length(line_input)*line_osr);
    for i =1 :length(line_input)
        if(i<length(line_input))
            add_step = (line_input(i+1) - line_input(i) ) / line_osr;
        end
        for j = 1 : line_osr
            if(j==1)
                line_y((i-1)*line_osr+j) = line_input(i);
            else
                line_y((i-1)*line_osr+j) = line_y((i-1)*line_osr+j-1) + add_step;
            end
        end
    end
    fft_plot_tittle=sprintf("线性插值滤波器输出数据,OSR=%d\n",fft_osr);
    plotCustomSpectrum(line_y(pn*fft_osr/4+1-fft_osr+1 : pn*fft_osr-fft_osr+1),length(line_y)*3/4,Fs*fft_osr,fft_plot_tittle);%% 对数据进行FFT之前进行了截取,以避免频谱泄露

%%
function [freq, spectrum] = plotCustomSpectrum(data, N,Fs,plot_title)
    % PLOTCUSTOMSPECTRUM 计算信号频谱并绘图
    % 输入：
    %   data      - 输入信号（一维向量）
    %   duration  - 波形时长（秒）
    % 输出：
    %   freq      - 频率向量（Hz）
    %   spectrum  - 频谱幅度

    % 计算采样频率
    % N = length(data);
    % Fs = N / duration; % 采样频率 = 点数 / 时长

    % FFT计算
%     Y = fft(data,N);
%     data = data .*blackman(N);
    Y = fft(data);


    % 幅度谱计算
    P2 = abs(Y / N);          % 双边幅度谱
    P1 = P2(1:floor(N/2)+1); % 单边幅度谱
    P1(2:end-1) = 2 * P1(2:end-1); % 能量校正[3,9](@ref)
    spectrum = P1;

    % 生成频率轴
    freq = Fs * (0:length(spectrum)-1) / N;

    % 转换为dB单位
    spectrum = 20 * log10(spectrum );

    % ==== 绘图 ====
    figure;
    hold on;

    % ==== 主瓣和噪声划分 ====
    mainlobe_mask = freq <= 20000;  % 主瓣定义为 0~20kHz
    noise_mask    = ~mainlobe_mask;

    % 主瓣用红色线绘制
    plot(freq(mainlobe_mask), spectrum(mainlobe_mask), 'r', 'LineWidth', 1.5);
    % 噪声用蓝色线绘制
    plot(freq(noise_mask), spectrum(noise_mask), 'b--', 'LineWidth', 1);
    % 画辅助线标注主瓣界限
    ylimits = ylim;
    line([22000 22000], ylimits, 'Color', [0.5 0.5 0.5], 'LineStyle', ':');

    % 图形设置
    xlabel('频率 (Hz)');
    ylabel('幅度 (dB)');
    title(sprintf('%s\n主瓣 vs 噪声对比（主瓣 ≤ 20KHz）', plot_title));
    legend('主瓣频段 (≤20KHz)', '噪声频段 (>20KHz)', 'Location', 'best');
    xlim([0, min(Fs/2, 100e3)]);
    grid on;
end
%%

% 叠加一个白噪声
%     m=0:M-1;
%     sigma=mean(sin_wave.*sin_wave);%噪声方差
%     sin_wave= sin_wave';
%     v=sqrt(sigma).*randn([M,1]);%噪声系数
%     x=sin_wave+v;%叠加噪声的音频信号